Draft version February 2, 2008 

Preprint typeset using LATpjX style cmulateapj v. 11/12/01 



BINARY AND MULTIPLE STAR FORMATION IN MAGNETIC CLOUDS: BAR GROWTH AND 

FRAGMENTATION 

FUMITAKA NAKAMURA 

Faculty of Education and Human Sciences, Niigata University, 8050 Ikarashi-2, Niigata 950-2181, Japan; 

fnakamur@ed.niigata-u.ac.jp 
AND 



Zhi-Yun Li 

Department of Astronomy, University of Virginia, P.O. Box 3818, Charlottesville, VA 22903; 

zl4h@virginia.edu 
Draft version February 2, 2008 

ABSTRACT 

In the standard scenario of isolated low-mass star formation, strongly magnetized molecular clouds 
are envisioned to condense gradually into dense cores, driven by ambipolar diffusion. Once the cores 
become magnetically supercritical, they collapse to form stars. Previous studies based on this scenario 
are limited to axisymmetric calculations leading to single supercritical core formation. The assumption 
of axisymmctry has precluded a detailed investigation of cloud fragmentation, generally thought to 
be a necessary step in the formation of binary and multiple stars. In a series of papers, we studied 
the non-axisymmetric evolution of initially magnetically subcritical clouds, using a two-dimensional 
magnetohydrodynamic code based on the physically motivated thin-disk approximation. We found that 
such clouds become unstable to non-axisymmetric perturbations after the supercritical cores are formed 
due to ambipolar diffusion. In this paper, we focus on the evolution of clouds perturbed by anm=2 
mode of a modest fractional amplitude of 5%, with an eye on binary and multiple star formation. 
We show that for a wide range of initial cloud parameters, the m — 2 mode grows nonlinearly into a 
bar during the isothermal collapse after the supercritical core formation. The instability is driven by 
the domination of the magnetically-diluted gravity over the combined thermal and magnetic pressure 
gradient in the supercritical cores. Such gravity-dominated cores can break up into fragments during or 
after the isothermal phase of cloud evolution. 

The outcome of fragmentation depends on the initial cloud conditions, such as the magnetic field 
strength, rotation rate, amount of cloud mass (relative to thermal Jeans mass) and mass distribution. 
It is classified into three different types: (1) separate core formation, in which the bar (m = 2) mode 
breaks up into two separate cores during the isothermal collapse, with a core separation of order 10 4 
AU, (2) bar fragmentation, in which the m = 2 mode evolves into a needle-like, opaque "first bar" (at a 
density n ;> 10 12 cm~ 3 ), which breaks up into multiple fragments with initial masses of order 1CT 2 M Q 
and separations of order 10 2 -10 3 AU, and (3) disk fragmentation, in which the bar growth remains slow 
during the isothermal collapse and the central region evolves into a rapidly rotating, opaque "first disk" , 
which breaks up into several self-gravitating blobs with separations less than the disk size (<; 10 2 AU). 
These three types of fragmentation loosely correspond to the empirical classification of embedded binary 
and multiple systems of Looney, Mundy, & Welch, based on millimeter dust continuum observations. 
The well-studied starless core, L1544, appears to belong to the bar fragmentation type. We expect it to 
produce a highly elongated, opaque bar at the center in the future, which should break up into fragments 
of initial masses in the substellar regime. 

Subject headings: binaries: formation — ISM: clouds — ISM: magnetic fields — MHD — stars: 
formation 



1. INTRODUCTION 

Most main-sequence stars are found in binary or low- 
order multiple-star systems (Duquennoy & Mayor 1991; 
Fischer & Marcy 1992). Recent observations of star- 
forming regions have revealed that the binary frequency of 
pre-main sequence stars is as high as or even higher than 
that of the main-sequence stars (Leinert et al. 1993; Ghez 
et al. 1997; Kohler et al. 2000; Mathieu et al. 2000). For 
example, in the well-studied nearby low-mass star forming 
Taurus-Auriga region, the binary frequency of pre-main 
sequence stars is a factor of two higher than that of the 
main-sequence stars. Other star-forming regions such as 
Ophiucus, Lupus, Chamaelcon, Scorpius-Centaurus, and 



Corona Australis show a similar excess of binary fraction. 
Observations of binary stars in young star clusters such 
as Orion and Pleiades indicate that the binary frequency 
of the pre-main sequence stars is comparable to that of 
the main-sequence stars (Prosser et al. 1994; Padgett et 
al. 1997; Bouvier et al. 1997). These observations suggest 
that binary formation is the primary mode of star forma- 
tion and that binaries should form prior to the pre-main 
sequence phase, most likely through fragmentation dur- 
ing the gravitational collapse of the parent dense molec- 
ular cloud cores (Mathieu 1994; Bodcnhcimcr et al. 2000; 
White & Ghez 2000). Recent detections of binaries and 
multiple systems at the earliest observable phases of star 
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formation (e.g., Class and Class I) support the fragmen- 
tation scenario (e.g., Wootten 1989; Fuller, Ladd, & Ho- 
dapp 1996; Looney et al. 1997, 2000; Terebey & Padgett 
1997). How the fragmentation occurs is a long-standing, 
unresolved problem. 

Many theoretical studies of star formation have concen- 
trated on the formation of single stars (Shu, Adams, & 
Lizano 1987). In the standard picture of single star for- 
mation, a strongly magnetized dense molecular cloud core 
is taken as the initial state of star formation. The dense 
core gradually contracts as the magnetic support weakens 
in the high density part due to the ambipolar diffusion. 
When the magnetic field becomes too weak, it is over- 
whelmed by self-gravity, and the cloud begins a dynamic 
contraction to form a star. The process of magnetic leak- 
age during the cloud contraction has been extensively in- 
vestigated both analytically (e.g., Mestcl & Spitzer 1956; 
Nakano & Tademaru 1972) and numerically (e.g., Nakano 
1979; Black & Scott 1982; Lizano & Shu 1989; Ciolek & 
Mouschovias 1993; Basu & Mouschovias 1994; Li 2001). 
Specific models based on this scenario appear able to ex- 
plain the observations of well-studied dense cores Bl and 
L1544 (Crutcher et al. 1994; Ciolek & Basu 2000; Li et al. 
2002). However, these studies have assumed axisymmetry 
of the model cloud, and are unable to address the issue 
of fragmentation. On the other hand, numerical hydrody- 
namic calculations of cloud fragmentation have neglected 
the effects of magnetic fields (e.g., Tohline 2002), which 
are potentially important or even critical. 

To quantify the role of magnetic fields on cloud frag- 
mentation, we have made a start in investigating the 
non-axisymmetric evolution and fragmentation of initially 
magnetically subcritical clouds, taking into account the 
effect of magnetic diffusion. Our ultimate goal is to eluci- 
date the fragmentation process in magnetically supported 
clouds, and extend the current standard scenario of low- 
mass star formation to binaries and multiple systems. 

Magnetically subcritical clouds are stable to both 
dynamical collapse and fragmentation (e.g., Tomisaka, 
Ikeuchi, & Nakamura 1988). Obviously, it is difficult 
to form binary or multiple stars by direct fragmentation 
of such clouds. Three-dimensional magnetohydrodynamic 
(MHD) simulations of magnetic cloud collapse have shown 
that a frozen-in magnetic field stifles cloud fragmentation 
(Dorfi 1982; Benz 1984; Phillips 1986a,b). In lightly ion- 
ized molecular clouds, however, ambipolar diffusion can 
initiate dynamical contraction in magnetically supercriti- 
cal cores by redistributing the magnetic flux in the cloud. 
The supercritical cores can be unstable to fragmentation if 
they have masses well above the Jeans limit (e.g., Mestel 
& Spitzer 1956; Galli et al. 2001). 

Recently, Boss (2000, 2002) carried out three- 
dimensional collapse calculations of initially magnetically 
supported clouds, taking into account several properties 
of magnetic fields approximately. He concluded that the 
cloud fragmentation is enhanced by magnetic fields be- 
cause the magnetic tension allows more time for the per- 
turbations to grow off center nonlinearly, which presents 
the formation of a single, dominant object at the center. 
Li (2001) studied ambipolar diffusion-driven collapse of 
magnetized cloud cores with one-dimensional axisymmet- 
ric calculations but treating the magnetic forces and am- 



bipolar diffusion more accurately. He also found a simi- 
lar beneficial effect of magnetic tension. He classified the 
evolution of magnetically subcritical axisymmetric clouds 
into two types, depending mainly on the initial cloud mass 
and initial density distribution. When the initial cloud 
has a relatively low mass and/or a centrally-condensed 
density distribution, it collapses to form a single super- 
critical core (core-forming cloud). On the other hand, 
when the initial cloud is more massive (containing several 
Jeans masses or more) and has a relatively flat density 
distribution, it collapses to form a ring after the central 
region becomes magnetically supercritical (ring-forming 
cloud). One-dimensional axisymmetric evolution of the 
core-forming clouds has also been extensively investigated 
by Ciolek & Mouschovias (1993) and Basu & Mouschovias 
(1994). Li & Nakamura (2002) examined the ring-forming 
clouds and showed that they can fragment into multiple 
cores in the presence of non-axisymmetric density fluctua- 
tions during the isothermal collapse phase after the super- 
critical core formation. This type of ring fragmentation 
produces fragments of relatively large initial separations 
(<~ 10 4 AU). It may be responsible for the formation of 
small stellar groups such as those in the Serpens cloud 
core (Williams & Myers 2000) and aggregates of pre-stellar 
cores such as those in the p Oph B2 core (Motte, Andre, 
& Neri 1998) and around IRAS 03256+3055 in Perseus 
(Young et al. 2002). These studies illustrate the crucial 
role of magnetic diffusion in the fragmentation of magnetic 
clouds, with strong implications for binary formation. 

In this paper, we focus on the non-axisymmetric evo- 
lution of the core-forming clouds, with an eye of the for- 
mation of binary stars. Initial results are presented in 
Nakamura & Li (2002), where we demonstrated that the 
core-forming clouds become unstable to a bar mode af- 
ter the supercritical core formation [see also Nakamura 
& Hanawa (1997) for initially magnetically supercritical 
clouds]. Here, we extend the calculations to a larger pa- 
rameter regime. First, we briefly describe our numerical 
method and model in §2. In §3, we examine the depen- 
dence of bar growth on the initial parameters such as mag- 
netic strength, initial cloud mass, density profile, and ro- 
tation. We show that the core-forming clouds become un- 
stable to a bar mode during the dynamic collapse phase 
over a wide range of initial model parameters. Subsequent 
fragmentation process is investigated in §4. Based on the 
numerical results, we classify the fragmentation process 
into three types, which agree with the empirical classifi- 
cation of young binary and multiple systems by Looney 
et al. (2000) based on high-resolution millimeter observa- 
tions. This general agreement and other implications of 
our results on binary formation are discussed in §5. 

2. FORMULATION OF THE PROBLEM 

A detailed description of the model cloud and the prob- 
lem formulation are presented in §2 of Li & Nakamura 
(2002). Here, we give a brief summary. We consider 
an initially axisymmetric, magnetically subcritical model 
cloud that is in a mechanical equilibrium, with a disk-like 
mass distribution due to settling along ordered magnetic 
field lines. In such a subcritical cloud, force balance along 
the field lines is achieved during the quasi-static phase of 
evolution and is well maintained even during the subse- 
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quent dynamic collapse (e.g., Fiedler & Mouschovias 1993; 
Nakamura, Hanawa, & Nakano 1995). Thus, we adopt 
the standard thin-disk approximation in which exact force 
balance along the vertical direction is assumed during the 
entire evolution (e.g., Ciolek & Mouschovias 1993; Basu & 
Mouschovias 1994; Shu & Li 1997; Nakamura & Hanawa 
1997; Stchlc & Spruit 2001; Li 2001). The magnetohy- 
drodynamical equations are integrated in the vertical di- 
rection and the cloud evolution is followed in the (x-y) 
disk plane of a Cartesian coordinate system (x, y, z) with 
a two-dimensional MHD code. The distribution of mag- 
netic field is computed in three-dimensional space under 
the assumption that the magnetic field is current-free out- 
side the disk. 

We assume an isothermal equation of state below a di- 
mensionless surface density S cr , and denote the isothermal 
sound speed by c s . Beyond S cr , we change the equation of 
state from isothermal to adiabatic, with an index of 5/3, 
to mimic the transition to optically thick regime of cloud 
evolution. In this paper, we set S cr = 1.9 x 10 4 £ ,rcf , cor- 
responding to the volume density 

Per ^ 10^ cm ^ (see eq. 
[1] below for the definition of S , ro f). 

The initial conditions for star formation are not well 
determined either observationally or theoretically. Follow- 
ing Basu & Mouschovias (1994) and Li (2001), we pre- 
scribe a slowly-rotating, axisymmetric reference state for 
our model cloud, with the distributions of mass, magnetic 
field, and velocity given by 

S ref (a:,y) = - S °' rcf , (1) 



[1 + (r/r )«] 4/n 



-B Ziro f(x,y) 
V<j,{x,y) = 
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= const. 
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(2) 
(3) 



r o + V r o + r 

where r = (x 2 + y 2 ) 1 ^ 2 is the distance from the origin, 
Cms = c,s(l +LQ) 1 / 2 is essentially the magnetosonic speed, 
and the parameter lo measures the rotation rate. The ex- 
ponent n controls the surface density profile and (indi- 
rectly) the amount of mass in the central plateau region 
where the distribution is more or less uniform. The back- 
ground field strength B ZyOQ is characterized by the dimen- 

sionless flux-to-mass ratio r = B Z;00 /(27rG 1,/2 £o,ref)- The 
reference clouds are not in an equilibrium state and are 
allowed to evolve into one with magnetic field frozen-in, 
before ambipolar diffusion is turned on at time t = 0. To 
save computational time, we obtain the equilibrium state 
using the one-dimensional axisymmetric MHD code of Li 
(2001), which gives results indistinguishable from those 
from the two-dimensional code assuming axisymmetry. 

At t = 0, we impose on top of the equilibrium surface 
density Y^ (x,y) anm = 2 perturbation of relative ampli- 
tude A, 

Z(x,y) = Z {x,y) [l + A cos(20)], (4) 

where (f> is the azimuthal angle measured from the x-axis. 
The subsequent, ambipolar diffusion-driven evolution of 
the perturbed cloud is followed numerically with the two 
dimensional MHD code. As a boundary condition, we 
fixed p and B z at the cloud outer radius, taken to be twice 
the characteristic radius r . 

The computations are carried out using non-dimensional 
quantities. The units we adopted are c s for speed, X , r ef 



for surface density, 27rGEo.rof for gravitational accelera- 
tion, and £?z. oo for magnetic field strength. The units for 
length and time are, respectively, L = c 2 / (27rG£ ,ref ) and 
to = c s /(27rGSo,rof)- They have typical values of 



Lo = 8.43 x 10 2 pc ( — I s )V , , 

V0.33 km s L / \10 1 g cm 



^0,rcf 



to = 2.50 x 10 5 yr f — — ^ T ) 

3 V0.33 kms" 1 / 



33 km s- 1 / V 10~ 2 g cm 



So.rof 



2 / ' 



(5) 
(6) 



where c s = 0.33 km s _1 corresponds to an effective tem- 
perature of T c s = 30 K. 

The initial state of our model cloud is completely spec- 
ified by four parameters: the characteristic radius r , ex- 
ponent n in the reference surface density distribution, di- 
mcnsionless flux-to-mass ratio To, and rotation parame- 
ter lo. We adopt a relative small rotation parameter of 
lo <J 0.2 — 0.3, consistent with the slow rates of rotation 
observed in molecular cloud cores (Goodman ct al. 1993). 
The slow rotation is thought to be due to magnetic braking 
by an ambient medium (Nakano & Tademaru 1972; Basu 
& Mouschovias 1994), which is not treated in our paper. 
We comment briefly on the potential effects of magnetic 
braking in § 5.3. 

Besides the four parameters characterizing the initial 
cloud, a parameter A is needed to specify the fractional 
amplitude of the m = 2 bar mode perturbation applied 
after the reference cloud settles into the equilibrium con- 
figuration. We have carried out a systematic parameter 
survey and selected the models listed in Table 1 to illus- 
trate the salient features of magnetic cloud evolution. The 
numerical results are described in the next two sections (§3 
and §4). 

3. BAR GROWTH IN MAGNETICALLY SUBCRITICAL 
CLOUDS 

As mentioned earlier, Li (2001) classified the evolution 
of magnetically subcritical axisymmetric clouds into two 
types, depending mainly on the initial cloud mass and 
initial density distribution: (a) core-forming cloud, which 
collapses to form a single supercritical core and (b) ring- 
forming cloud in which a supercritical ring is formed after 
the central region becomes supercritical owing to ambipo- 
lar diffusion. The ring fragments into multiple supercrit- 
ical cores in the presence of density fluctuations during 
the isothermal collapse phase (Li & Nakamura 2002). In 
this paper, we focus on the core-forming clouds and show 
that these clouds become unstable to bar formation af- 
ter the supercritical core formation over a wide range of 
initial model parameters. The subsequent evolution and 
fragmentation of the bar are discussed in §4. 

3.1. Non-Rotating Cloud 

To begin with, we consider a non-rotating (lo = 0) cloud 
with the reference state specified by the set of parameters 
r = 10tt, n = 2, and r = 1.5. The cloud is therefore 
magnetically subcritical, with the field strength 50% above 
the critical value at the center. The cloud so specified pro- 
duces a dense supercritical core in the absence of any non- 
axisymmetric perturbations. The reference cloud is al- 
lowed to evolve into an equilibrium configuration, with the 
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magnetic field frozen-in. After the equilibrium is reached, 
we reset the time to t = 0, turn on ambipolar diffusion, 
and add an m = 2 perturbation to the surface density dis- 
tribution, with a fractional amplitude of A = 0.05. The 
evolution of such a non-axisymmctrically perturbed cloud 
is followed numerically to progressively higher densities 
and smaller scales. The results are shown in Fig. 1 and 
described as follows. 

The evolution of the model is qualitatively similar to 
that of the model shown in Fig. 1 of Nakamura & Li 
(2002), where a different reference density distribution 
(less centrally-condensed than the model shown in this 
subsection) is adopted. The evolution is followed to a 
higher density here. As is usual with ambipolar diffusion- 
driven evolution, the cloud spends most of its time in the 
subcritical phase when the dimensionlcss minimum flux- 
to-mass ratio (at the density peak) is greater than unity, 
r c > 1 (see panels [a] and [b] of Fig. 1). During this pe- 
riod, gas motions are very subsonic and the central part 
of the cloud oscillates with a small amplitude. The os- 
cillation indicates that the cloud is stable to the m = 2 
perturbation during the subcritical phase, in agreement 
with linear analysis by Nakano (1988) and others. Indeed, 
we have verified numerically that for a frozen-in field the 
subcritical cloud remains stable to this and other higher 
(e.g., m = 3, 4, . . .) modes. (We followed the evolution of 
the frozen-in model for several periods of oscillation and 
confirmed that the cloud oscillates with almost the same 
amplitude as the initial.) Once the minimum flux-to-mass 
ratio has dropped below the critical value, the contrac- 
tion becomes dynamic and the infall speed eventually ex- 
ceeds the sound speed at the center (see panels [c] and 
[d]). By the time shown in panel (c), the m = 2 mode 
has grown significantly, resulting in a bar-like core at the 
center with an aspect ratio of roughly 2. The bar growth 
remains relatively modest (panel [d]) until the very end 
of the runaway collapse, when the growth rate increases 
dramatically (panel [e]). Beyond a dimensionlcss critical 
surface density of 1.9 x 10 4 , corresponding to a volume 
density of n# <~ 10 12 cm -3 , we change the equation of 
state from isothermal to adiabatic, with an index of 5/3, 
to mimic the transition to optically thick regime of cloud 
evolution. When the maximum surface density reaches 
~ 10 5 , the contraction along the minor axis of the bar de- 
celerates significantly and an accretion shock is formed at 
the surface of the bar. This shock-bound bar is analogous 
to the "first core" of spherical calculations, and will be re- 
ferred to as the "first bar" below. Immediately after the 
shock formation, the aspect ratio of the first bar contin- 
ues to increase. It reaches a maximum of about 10 by the 
time shown in panel (e). Thereafter, the infall along the 
major axis becomes significant, forming a round core at 
the center. The aspect ratio of the central round core ap- 
proaches unity as the collapse proceeds. Since the round 
core resembles the first core of the spherical calculations, 
we call it the "truly first core". As will be shown in §4, 
this round core evolves into a rapidly-rotating disk in the 
presence of an initial, slow cloud rotation. 

To quantify the geometrical change of the supercritical 
core, we plot in Figure 2 the aspect ratio of the darkest 
region in Figure 1 (with £ > 10 _1 / 2 S max ) as a function 
of the maximum surface density, S max , together with the 



flux-to-mass-ratio T at the center. The aspect ratio is de- 
fined as R = a y /a x where a x and a y are, respectively, the 
effective scale heights along the x and y axes, and evalu- 
ated as 

a x = (m~ x [ Y,x 2 dxdy) 1 , (7) 

V t /£>10- 1 /2£ max / 

and 

a v = ( M" 1 f Xy 2 dxdy) , (8) 

V ■/E>10- 1 /2E max / 

where M is the total mass in the region with £ > 
10 _1 / 2 S max . During the subcritical phase, the aspect ra- 
tio of the bar does not grow in time; it oscillates around 
unity, signifying that the subcritical cloud is stable to the 
m — 2 mode. As the central region becomes more and 
more supercritical, the bar mode begins to grow signifi- 
cantly. During this supercritical period, the aspect ratio 
of the bar stays more or less frozen at R <~ 1.5 before the 
central density reaches £ ~ 10 3 . This freezing of aspect 
ratio is related to the dynamical state of the disk. As 
shown in §3.4, the mass of the central flat region is merely 
~ 1.5 times the effective Jeans mass (as measured by the 
central surface density, including the magnetic contribu- 
tion). Such a relatively low-mass disk is only marginally 
unstable to the non-axisymmetric perturbation. In other 
words, the time scale of bar growth is comparable to that 
of dynamic contraction. The flux-to-mass ratio does not 
change much either during this period, as a result of the 
rapid collapse which prevents the magnetic flux from leak- 
ing out significantly. The bar growth speeds up toward the 
very end of the supercritical collapse when the mass of the 
central flat region exceeds about twice the effective Jeans 
mass. It continues after the first bar formation. When 
the central surface density reaches 10 5 , the infall motion 
along the major axis of the bar becomes significant and 
the aspect ratio begins to drop steeply to unity, forming 
the truly first core. 

To study the structure of the bar, we show in Figure 
3 the distributions of the surface density, magnetic field, 
and infall velocity along the axes as a function of the dis- 
tance from the density maximum, at seven different times. 
The solid and dashed curves indicate the distributions of 
the physical quantities along the major (y) and minor (x) 
axes, respectively. During the early phase of supercriti- 
cal collapse when the aspect ratio of the bar is frozen at 
R ~ 1.5, the surface density distributions along both axes 
qualitatively resemble the self-similar distribution of the 
axisymmetric collapse; they are flat near the density max- 
imum and can be approximated by a power-law of £ cx dr 1 
at the magnetically supercritical envelope, where d is the 
distance from the density peak. The X oc c? _1 distribution 
corresponds to a volume density distribution of p oc d~ 2 
because of the force balance in the vertical direction (which 
yields p oc £ 2 ). While the central density increases rapidly, 
the surface density distributions change little in the enve- 
lope. The central flat region shrinks in size in accordance 
with the increase in the central density. The distributions 
of the vertical field B z are qualitatively similar to those 
of £ in the magnetically supercritical region and almost 
proportional to the local surface density (B z oc £). The 
x- and y-components of the magnetic field (in the disk 
plane) increase approximately linearly with distance d in 
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the central flat region, peaking near the edge of the flat 
region, before falling off as B x oc d^ 1 and B y oc d^ 1 . In 
the magnetically supercritical envelope, B x and B y are 
comparable to B z , indicating that the field lines near the 
disk plane incline with an angle of ~ 45 degree with re- 
spect to the disk plane. The angle approaches ~ 90 degree 
in the outer, magnetically subcritical envelope. The infall 
speeds are proportional to the distance d near the center, 
becoming more or less flat in the supercritical part of the 
envelope, before falling to subsonic values in the outer sub- 
critical envelope. At S max = 10 4 , the infall speeds in the 
supercritical envelope reach ~ 2c s . 

As the collapse proceeds and the bar aspect ratio be- 
comes large, the distributions of physical quantities grad- 
ually deviate from the self-similar distributions. For ex- 
ample, the surface density and vertical magnetic field 
(B z ) profiles along the major axis become steeper as the 
bar elongates rapidly. When the central surface density 
reaches ~ 2x 10 4 , a shock is formed at x ~ 10~ 4 ro. There- 
after, the infall speed rises in the envelope of the bar. The 
infall speed along the minor axis follows v x oc d -1 / 2 just 
outside the first bar, indicating a free-fall collapse towards 
the bar. 

Interestingly, after the bar becomes opaque to dust emis- 
sion, a local density peak is formed at each end of the first 
bar. This density peak formation is due to the "edge" 
effect, which was first calculated by Larson (1972) and 
then extensively investigated by Bastien (1983) and Bon- 
nell & Bastien (1991). In this model, the density peaks do 
not grow nonlinearly; they tend to disappear in the sub- 
sequent collapse because the central truly first core grows 
more rapidly in mass owing to the supersonic infall along 
the major axis, which dominates the gravitational poten- 
tial in the high density region. In §4, we show that for 
more massive clouds in which the aspect ratios of the su- 
percritical cores are larger, the density peaks can grow in 
time because the infall speeds along the major axis remain 
subsonic. 

The model discussed in this subsection serves as a stan- 
dard against which other models will be compared. 

3.2. Dependence of Model Parameters: T , r , and n 

The properties of the bar growth depend on the ini- 
tial model parameters. Here, we investigate the effects of 
the initial flux-to-mass ratio To, the characteristic radius 
ro, and the density profile parameter n for non-rotating 
clouds. The effect of rotation is examined in separate sub- 
sections §3.3 and §4.1. Exact values of T are difficult to 
determine observationally; they probably lie within a fac- 
tor of two of unity (Crutcher 1999), after correcting for 
likely projection effects (Shu et al. 1999). Thus, we inves- 
tigated the dependence on r in the range of 1 <; r <; 3. 
The characteristic radius ro is to be compared with 2n, the 
(dimcnsionlcss) critical wavelength for the (thermal) Jeans 
instability in a disk of uniform mass distribution (Larson 
1985). The parameter n controls the profile of the cloud 
surface density distribution; those with smaller values of 
n are more centrally peaked. 

To examine the effect of the initial flux-to-mass ratio 
To, we first hold the cloud radius and profile parameter 

1 We note that the magnetic field has an opposite, beneficial effect on the fragmentation of ring-forming clouds. More strongly magnetized 
clouds contract more slowly, allowing more time for non-axisymmetric perturbations to grow off-center. As a result, multiple fragmentation is 
achieved more readily, at a lower surface density. See Li & Nakamura (2002) for details. 



fixed at the standard values ro = 107r and n = 2, and re- 
duce the flux-to-mass ratio r . Representative results are 
shown in panel (a) of Fig. 4, at a time when the maxi- 
mum surface density E max = 10 4 . Panel (a) is the snap- 
shot of a cloud with Tq — 1.25 (a weaker initial magnetic 
field), perturbed by an m = 2 mode of the standard frac- 
tional amplitude A = 0.05. Comparison with panel (d) of 
Fig. 1 where the standard case at the same S max is plot- 
ted shows that a weaker magnetic field promotes the bar 
growth. This slower bar growth in the stronger field case 
is quantitatively shown in panel (a) of Figure 2 (compare 
the solid curve with the dotted-dashed curve). The reason 
for this behavior appears to be that, for the cloud with 
a stronger initial magnetic field, the magnetic support re- 
mains more important even in the supercritical phase. In 
other words, the averaged flux-to-mass ratio in the dark- 
est region (S > 10 _1 / 2 X max ) is larger for the stronger field 
model. The stronger field tends to retard the deformation 
of the supercritical core. We have confirmed this negative 
effect of the magnetic field on bar growth by following the 
evolution of the clouds with r = 2 and 3. 1 

We next consider the effect of the characteristic radius 
ro which controls the total cloud mass. Increasing the ra- 
dius r is expected to have a positive effect on bar growth. 
This is indeed the illustrated in panel (b) of Fig. 

4, where the snapshot of a cloud with r = 157T, 1.5 times 
the standard value, is shown at S max = 10 4 . The aspect 
ratio of the bar reaches ~ 7 at the time when the surface 
density reaches 10 4 . 

The profile parameter n has an effect similar to the 
parameter ro. For instance, a reduction in n would de- 
crease the number of Jeans masses contained in the cen- 
tral plateau region. We showed in a previous paper (Li & 
Nakamura 2002) that for the ring-forming case, the clouds 
with more peaked density profiles are less prone to frag- 
mentation. This negative effect of the centrally condensed 
density distribution is also observed for non-magnetized 
clouds (Boss 1993). We find the same trend for the core- 
forming clouds as well. The aspect ratio of the bar in- 
creases more rapidly for larger n. This is illustrated in 
panel (c) of Fig. 4, where n = 4, twice the standard value, 
is chosen. The bar is more elongated than the standard 
case, even though the characteristic radius r = 5tt is only 
half of the standard value. 

3.3. Slowly-Rotating Cloud 

To gauge the effect of slow rotation, we repeat the evo- 
lution of a cloud with all parameters identical to those of 
the standard model in §3.1 except the rotation parameter 
u>, which is now set to 0.1 instead of zero. Snapshots of 
the surface density distribution and velocity field at three 
selected times are shown in Fig. 5. The general trend of 
the cloud evolution is similar to that in the non-rotating 
case; in the subcritical phase, the cloud is stable to non- 
axisymmetric perturbations and the m = 2 mode begins 
to grow in time only after the central region becomes mag- 
netically supercritical. Comparison between Fig. Id and 
Fig. 5b indicates a beneficial effect of the slow rotation on 
the bar growth. More quantitatively, the aspect ratio of 



6 



Nakamura & Li 



the bar is always larger for the rotating cloud at the time 
when the maximum density reaches the same value, before 
the bar becomes opaque to dust emission. This effect de- 
rives from the fact that rotation retards the radial infall, 
which allows more time for perturbations to grow. 

In §4.1 below, we will show that the slow-down effect 
due to rotation plays a dramatic role in the cloud fragmen- 
tation for more rapidly-rotating, less centrally-condensed 
clouds, which can fragment into two cores during the 
isothermal collapse phase soon after the supercritical core 
formation. 

3.4. Dynamical State of Supercritical Cores: 
Axisymmetric Collapse 

The bar growth is closely related to the dynamical state 
of the supercritical core. In this subsection, to clarify the 
physics of the bar growth, we follow the evolution of un- 
perturbed axisymmetric clouds with the one-dimensional 
code and measure the dynamical forces in the supercriti- 
cal cores. For simplicity, we do not change the equation 
of state for the one-dimensional calculations even when 
the surface density exceeds the critical value of 1.9 x 10 4 . 
Previous axisymmetric calculations of ambipolar diffusion- 
driven cloud evolution have shown that the collapse of su- 
percritical cores proceeds self-similarly: the distributions 
of the surface density and magnetic field are almost con- 
stant in the central flat region and decrease with radius 
as £ oc r -1 and B z oc r _1 (Ciolek & Mouschovias 1993; 
Basu & Mouschovias 1994). In the high density region 
of the supercritical core, the decrease of the flux-to-mass 
ratio levels off (at T ~ 0.4) due to the rapid dynamic 
collapse which prevents significant magnetic flux leakage. 
This is shown in panels (a) and (c) of Figure 6, where the 
flux-to-mass ratio normalized by the critical value is de- 
picted as a function of radius for two representative mod- 
els. Shu & Li (1997) and Li & Shu (1997) called such a 
cloud with a constant flux-to-mass ratio "isopedic" . In the 
isopedic cloud, the gravity (F grav = T,g r ) and gas pressure 
(-Fgas = c 2 s dY>/dr) are proportional to the magnetic tension 
(F ten = B r B z /2ir) and pressure (F B = H d(B 2 / An) / dr) , 
respectively (Shu & Li 1997; Li & Shu 1997; Nakamura & 
Hanawa 1997), where H is the half thickness of the disk 
(see eq. [5] of Li & Nakamura (2002) for the definition of 
H). In fact, the supercritical cores have such a character- 
istic. 

Figure 6 shows the evolution of the ratio of the effective 
gravity to the effective pressure at six different times for 
two models. Here, the effective pressure and gravity are 
defined as F gas + and F grav + F ton , respectively. At the 
initial state (t = 0) , this ratio is equal to unity, indicating a 
force balance in the initial cloud. As the collapse proceeds, 
this ratio increases due to ambipolar diffusion, which re- 
duces the magnetic support in the central high density 
region. When the central surface density is in the range of 
10 3 — 10 4 , the effective gravity is larger than the effective 
pressure, but only by a factor of about 1.5 — 2. This cor- 
responds to a mass of the central flat region only about 1.5 
— 2 times the effective Jeans mass. Such a relatively low- 
mass core is only marginally unstable to non-axisymmetric 
perturbations, which is probably the reason why the as- 
pect ratio of the bar remains almost frozen at a modest 
value of 1.5 — 2 in the presence of an m = 2 perturba- 



tion. When the force ratio exceeds ~ 2, the bar growth is 
accelerated remarkably by the Lin-Mestel-Shu type grav- 
itational instability (Lin, Mestel, & Shu 1965), as seen in 
Fig. 2. For both clouds, the force ratios approach a higher 
value of 3 — 4 as the collapse proceeds, although a conver- 
gence hasn't been reached by the time the central surface 
density reaches the critical value of 1.9 x 10 4 . (The con- 
vergence toward a self-similar solution is slower for more 
centrally-condensed clouds.) This behavior suggests that 
the isothermal evolution tends toward a self-similar solu- 
tion. 

A self-similar solution reproducing the pre-protostellar 
collapse of an isopedic magnetic disk is derived approxi- 
mately by Nakamura & Hanawa (1997), who considered a 
frozen-in magnetic field [see also Saigo & Hanawa (1998) 
for the nonmagnctized case]. Basu (1997) constructed 
a self-similar collapse model including the effect of am- 
bipolar diffusion [see Li (1998) for the self-similar collapse 
after the protostellar core formation]. The behavior of 
his solution is in good agreement with that of Nakamura 
& Hanawa (1997) in the late, self-similar phase of dy- 
namic collapse, when the effective gravity is about 3 — 4 
times the effective pressure force and the asymptotic ve- 
locity is about twice the effective sound speed. Naka- 
mura & Hanawa (1997) investigated the dynamical stabil- 
ity of the self-similar solution to non-axisymmetric pertur- 
bations and found that the self-similar solution is unstable 
to the to = 2 mode and stable to the higher modes. This 
characteristic is consistent with the numerical results pre- 
sented in the previous subsections. Thus, we conclude that 
the tendency for the supercritical collapse to approach the 
gravity-dominated self-similar solution is responsible for 
the bar formation during the dynamic collapse. 

4. FRAGMENTATION OF SUPERCRITICAL CORES 

In §3, we demonstrated that initially magnetically sub- 
critical clouds tend to be unstable to an to = 2 mode, 
forming bar-like cores. In this section, we investigate the 
subsequent evolution of the bar mode and show that the 
elongated supercritical cores can break up into discrete 
fragments. We classify the fragmentation process into 
three types, depending on the initial model parameters: 
(1) separate core formation, in which the to = 2 mode 
breaks up into two separate cores in the isothermal col- 
lapse phase, (2) bar fragmentation, in which the first bar 
breaks up into multiple fragments that move mainly along 
the major axis of the bar, with frequent merging among the 
fragments, and (3) disk fragmentation, in which the central 
region evolves into a rapidly-rotating disk where another 
bar mode is excited, fragmenting into blobs. In the follow- 
ing, we select three models to illustrate the salient features 
of these cases. 

4.1. Separate Core Formation: Effect of Faster Rotation 

In §3.3, we showed that rotation slows down the dy- 
namic collapse, allowing more time for perturbations to 
grow. Here, we demonstrate that a relatively fast rota- 
tion can affect the outcome of cloud evolution dramati- 
cally, particularly for n J> 4 clouds whose initial parameter 
range is intermediate between those for core-forming and 
ring-forming clouds in the axisymmetric calculations (Li 
2001). Panels (b) and (c) of Fig. 7 show the snapshots of 
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the surface density distribution and velocity field at two 
selected times for the model with n = 4, u = 0.2, r = 57r, 
and T = 1.5. This cloud is categorized as a core- forming 
cloud from the axisymmetric one-dimensional calculation. 
A cloud larger in radius by only a few tens of percent would 
collapse to form a ring instead of a single core. For com- 
parison, the snapshot of a more slowly rotating cloud with 
cj = 0.1 is shown in panel (a) at the time when S max = 10 2 . 
The model with u> = 0.1 shows the evolution qualitatively 
similar to the model in §3.3, i.e., after the supercritical 
core formation, a bar mode grows in time without frag- 
mentation. On the other hand, when the rotation rate is 
doubled (ui = 0.2), the bar grows more rapidly even at the 
early time shown in panel (b) , when the maximum surface 
density reaches 10. By this time, the entire bar (the dark- 
est part in panel [b]) becomes magnetically supercritical. 
By the time shown in panel (c) (S max = 10 2 ), the m = 2 
mode breaks up into two separate cores, each of which 
continues to collapse dynamically. (We checked that this 
model is stable to higher to modes.) The subsequent evo- 
lution of each core is qualitatively similar to that of the 
models shown in the previous section, i.e., the aspect ratio 
of each core is more or less frozen at ^ 2 during the early 
phase of supercritical collapse, and increases rapidly at the 
very end of the supercritical collapse. 

For the separate core formation case, the magnetic field 
has a positive effect on the fragmentation. We calculated 
the models of a weaker and a stronger magnetic field for 
the same set of parameters n — 4, ro = 5ir, lu = 0.2, and 
A = 0.05, and found that the cloud with a weaker magnetic 
field does not break up into two cores (it collapses into 
a single bar instead), whereas the cloud with a stronger 
magnetic field fragments into two cores at a lower density 
and with a wider separation. More centrally-condensed 
clouds with n = 2 and r <; 207r do not fragment into two 
cores even in the presence of a rapid rotation (w ~ 0.3). 
These dependences on r and n are the same as those 
of the multiple fragmentation case investigated by Li & 
Nakamura (2002). In this sense, the separate core case is 
similar to the multiple fragmentation case, although the 
bar fragmentation into two cores appears to depend more 
strongly on the rotation rate. 

At the time shown in panel (c), the separation of two 
cores is about a tenth of the initial cloud radius, which cor- 
responds to 5 x 10 4 AU with the fiducial choice of So = 0.01 
g cm' 2 and T e g = 30 K. This separate core formation case 
may thus be responsible for the formation of wide binaries 
with separations of the order of 10 4 AU, although substan- 
tial orbital evolution is likely after the initial fragmenta- 
tion, since the cores are not kept apart by rotation. 

To summarize, when the initial rotation is relatively fast 
and the density distribution in the central region is rela- 
tively flat, the to = 2 mode can break up into two sepa- 
rate cores in the isothermal collapse phase after the cen- 
tral region becomes supercritical (e.g., for n = 4, u> ~ 
0.2, T ~ 1.5, and r J> 5n). We note that this separate 
core formation differs from that found previously for non- 
magnetic clouds, which requires a rapid rotation or highly 
unstable initial condition (e.g., Tsuribe & Inutsuka 1999). 
In our case, the cloud is initially stable to gravitational 
collapse, and the fragmentation is induced by ambipolar 
diffusion. Indeed, we find that clouds with relatively weak 



magnetic fields do not fragment into separate cores during 
the isothermal collapse phase. 

4.2. Bar Fragmentation and Merging of Fragments 

According to linear perturbation theory, long filamen- 
tary clouds tend to break up gravitationally into fragments 
if their lengths exceed several times the diameter. This is 
indeed the case for the first bars formed in relatively mas- 
sive, centrally-condensed clouds (with a relatively large 
r and n <; 4, e.g., for n = 2, r ~ 1.5, u ~ 0.2, and 
ro ;> 157r), as illustrated in Fig. 8, where the snapshots of 
a cloud with n = 2, Tq = 1.5, ro = 157r, and u> = 0.2 are 
shown at three different times. In this model, fragmen- 
tation starts when the bar becomes opaque to dust emis- 
sion, which retards the radial contraction (at S c J> 10 5 ). 
At first, the bar produces two dense fragments near the 
center at £ ~ 10 5 , and then the fragmentation apparently 
propagates towards the ends of the bar. This apparent 
propagation of the fragmentation is due to the fact that 
the growth rate of the gravitational fragmentation is larger 
for the higher density region closer to the center. By the 
time shown in panel (a), there is a fragment formed at each 
of the two ends of the bar. This density peak formation 
at the end was first noticed by Larson (1972) and sub- 
sequently investigated extensively by Bastien (1983) and 
Bonncll & Bastien (1991). It is caused by a local minimum 
in the gravitational potential. The masses of the fragments 
formed by bar fragmentation are on the order of 1O _2 M0. 
These masses agree with those expected from the linear 
theory of fragmentation of cylindrical clouds (e.g., Larson 
1985; Nakamura, Hanawa, & Nakano 1993). By the time 
shown in panel (b), the two pairs of fragments near the 
center have merged to produce one, more massive, pair. 
This process of successive bar fragmentation and subse- 
quent merging is repeated, as shown in panel (c), where 
the two most massive fragments are about to merge at 
the center. We will discuss further the evolution of this 
multiple fragment system in §5. 

We note that bar fragmentation was also observed in the 
numerical simulations of non-magnetic cloud collapse (e.g., 
Boss 1993; Bate & Burkert 1997). The non-magnetic bars 
are typically much shorter than the one presented in this 
subsection. To form such a long bar in the non-magnetic 
case, the initial cloud must contain many thermal Jeans 
masses and thus be highly unstable. In our case, the ini- 
tial cloud is stablized by a strong magnetic field, whose 
support weakens with time as a result of ambipolar diffu- 
sion. Nakamura & Hanawa (1997) studied bar formation 
in clouds with frozen- in magnetic fields. They found bars 
that are prone to gravitational breakup according to the 
linear theory of fragmentation. These bars are however 
shorter than the one discussed here. They are expected to 
produce a smaller number of fragments. 

4.3. Disk Formation and Fragmentation 

When the initial clouds are centrally condensed and/or 
less massive (e.g., for n = 2, u <~ 0.1, r ~ 1.5, and 
ro <J 157t), the length of the bar does not exceed the criti- 
cal wavelength for gravitational fragmentation by the time 
of first bar formation. It would be difficult for such a short 
bar to break up into blobs gravitationally. In this case, a 
rapidly-rotating disk is formed in the presence of slow ro- 



8 



Nakamura & Li 



tation. Here, we show an example of the rotating disk 
formation. Figure 9 shows the snapshots at six different 
times of a model with n — 2, T = 1.5, r = 5ir, and 
lo = 0.1, which is identical to the slowly rotating model 
discussed in §3.3, except for a smaller characteristic ra- 
dius r (57r instead of 107r). For comparison, the linear 
density contours of the central disk is depicted in Fig. 10 
at three selected times. In this model, the bar mode does 
not grow much by the end of isothermal collapse because 
of the relatively small cloud radius tq 2 . Because the length 
of the bar is smaller than the critical wavelength for grav- 
itational fragmentation, it does not fragment. Instead, a 
round dense region, which we called the "truly first core" , 
is formed at the center due to a supersonic infall motion 
along the bar. 

Initially, the truly first core is supported primarily by 
the thermal pressure. As ambient gas with a larger an- 
gular momentum accretes onto it, the rotational support 
becomes dominant. As a result, the central core evolves 
into a rapidly-rotating, opaque disk, which we will some- 
times refer to as the "first disk" . The interaction between 
the disk and the large-scale bar mode outside the disk ex- 
cites another bar- mode growing inside the disk. We believe 
that the secondary bar growth is due to rotational insta- 
bility because the ratio of the rotational energy and the 
gravitational energy of the first disk reaches the critical 
value of f3 CI ~ 0.3 (e.g., Durisen et al. 1986; Bate 1998; 
Tohline 2002). By the time shown in panel (c) of Fig. 10, 
the bar breaks up into two blobs, which appear to be self- 
gravitating. We stopped our calculation at this epoch of 
fragmentation because the volume densities of the blobs 
reach the critical density beyond which H2 dissociation 
starts and the equation of state of the gas is expected to 
change again. 

The fate of the two blobs formed in the disk is uncer- 
tain. They may merge together during the subsequent 
evolution. Since the rotation period of the disk is faster 
than that of the large-scale bar outside the disk, the inter- 
action between the two bar modes generates several blobs 
in the disk, some of which appear to be self-gravitating. 
These blobs are expected to interact gravitationally and 
whether any of them can survive the interaction remains 
to be determined. This type of fragmentation and merging 
in a rotating disk was also observed in the numerical simu- 
lations of non-magnetic cloud collapse (e.g., Bonnell 1994; 
Whitworth et al. 1995; Burkert et al. 1997; Matsumoto & 
Hanawa 2002). 



suggest that strong magnetic fields tend to stifle fragmen- 
tation of dynamically collapsing clouds (Dorfi 1982; Bcnz 
1984; Philiips 1986a,b). For example, Phillips (1986a) 
performed three-dimensional MHD simulations on initially 
magnetically supercritical clouds and found that none of 
his clouds show any evidence for fragmentation even in the 
presence of a large density perturbation with a fractional 
amplitude of 50%. We emphasize that there are two im- 
portant differences between our model and the previous 
ones. First, previous studies considered a frozen-in field, 
whereas in our study ambipolar diffusion is a crucial in- 
gredient. Our initial cloud, which is supported by a strong 
magnetic field, is stable to gravitational fragmentation as 
well as dynamic collapse. The weakening of the magnetic 
support driven by ambipolar diffusion is what initiates the 
dynamic collapse and fragmentation. Second, the contri- 
bution of the magnetic tension to the cloud support is 
much weaker in the earlier studies; indeed, most of them 
took a uniform magnetic field as the initial state (which of 
course has no magnetic tension). The beneficial effect of 
magnetic tension on cloud fragmentation was discussed by 
Boss (2000, 2002), who showed that a strong magnetic ten- 
sion can slow down cloud contraction appreciably, which 
tends to enhance cloud fragmentation. 

Recently, Tohline (2002) reviewed the theoretical stud- 
ies on binary star formation, concentrating on the non- 
magnetic case. He found broad agreement among nu- 
merical simulations by various groups that non-magnetic 
clouds are susceptible to prompt fragmentation if they col- 
lapse from an initial configuration that is relatively uni- 
form in density and contains more than a few Jeans masses 
of material. If supported only by thermal pressure, such 
a configuration is highly unstable. The question is then 
(Tohline 2002): how is nature able to produce such an un- 
stable configuration in the first place? We believe that the 
key to resolve this apparent difficulty is the presence of 
a strong magnetic field, which stabilizes the entire cloud 
even when the cloud mass is well above one thermal Jeans 
mass 3 (e.g., Lizano & Shu 1989). Specific models of core 
formation in magnetically subcritical clouds appear able 
to explain the mass distribution and extended subsonic 
infall motions inferred in several starless cores (e.g., Lee, 
Myers & Taffala 2001; Ciolek & Basu 2000; Indebetouw & 
Zwcibel 2000) and are consistent with pattern of chemical 
differentiation observed in L1544 (Li et al. 2002), provided 
that the initial flux-to- mass ratio lies within a factor of two 
or so of the critical value, as adopted in this paper. 



5. DISCUSSION 

5.1. Fragmentation of Magnetic Clouds 

In the previous sections, we have shown that core- 
forming clouds tend to become unstable to the bar mode 
after the supercritical core formation. The numerical re- 
sults presented in this paper and in Li & Nakamura (2002) 
indicate that fragmentation of magnetic clouds happens 
readily over a wide range of initial cloud parameters. On 
the surface, this result appears to contradict earlier nu- 
merical studies of magnetic cloud fragmentation, which 

2 We estimate that the central plateau region of this cloud has a size roughly 2 — 3 times the effective Jeans length, judging from the fact that, 
for n = 2, reference clouds with tq < 2ir expand rather than contract. 

3 The cloud can also be supported by turbulence, although the details of cloud evolution driven by turbulence decay remain to be worked out 
(Myers 1999). 



5.2. Elongation of Dense Cloud Cores 

Dense cores of molecular clouds are observed to have 
significant elongation, with a typical aspect ratio of 2 in 
the plane of sky (Myers et al. 1991). Their true three- 
dimensional shapes are probably triaxial, as inferred sta- 
tistically by Jones, Basu & Dubinski (2001) and Goodwin, 
Ward- Thompson, & Whitworth (2002). We have previ- 
ously noted that the triaxial core shape can naturally re- 
sult from nonlinear growth of the bar mode in the direction 
perpendicular to the field lines, coupled with the flattening 
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of the mass distribution along the field lines (Nakamura & 
Li 2002). The models presented in §3 strengthen this con- 
clusion. We find that the aspect ratio of the bar is more 
or less independent of the initial model parameters. It 
stays at about 1.5 — 2 during the period when T decreases 
from ~ 1 to ~ 0.5, after the core has become supercritical 
(which makes bar growth possible) but before the ratio 
of effective gravity to effective pressure force exceeds ~ 2 
(when the bar elongation is greatly amplified through the 
Lin-Mcstel-Shu instability). This initial period of super- 
critical phase is also the phase most likely observed as a 
starless core, when the core is dense enough to be identified 
as a core using high density tracers (such as NH3 or CS), 
but before rapid collapse sets in which quickly produces a 
stellar object at the center, ending the starless phase. 

In this supercritical collapse phase, the decrease of the 
central flux-to-mass ratio saturates at about half the crit- 
ical value, almost independent of the initial parameters 
(see panel [b] of Fig. 2). This saturation of flux-to-mass 
ratio is also observed in axisymmetric models (Basu & 
Mouschovias 1995). It means that the supercritical cores 
remain significantly magnetized, and magnetic fields may 
not be completely ignored even during the dynamic col- 
lapse phase of star formation. 

5.3. Three Types of Binary and Multiple Star Formation 

Observations of nearby star forming regions have estab- 
lished that binary fraction among pre-main sequence stars 
is at least as high as among the main-sequence stars. In 
some nearby low-mass star forming regions like Taurus- 
Auriga, the binary frequency is even higher than that of 
main-sequence stars. In such low-mass star forming re- 
gions, stellar densities (n* ~ 10 pc~ 3 ) are much too low 
to allow the formation of binaries by capture or their de- 
struction by gravitational encounters within their ages of 
~ 1 Myr. This indicates that the observed high binary 
frequencies more or less reflect the pristine population of 
proto-binaries at the epoch of formation, suggesting that 
binary and small multiple stars are formed prior to pre- 
main sequence phase by gravitational fragmentation of 
parent dense molecular cloud cores. In fact, recent high 
resolution millimeter observations of dust continuum have 
revealed evidence of gravitational fragmentation of proto- 
stellar cores in the earliest, Class and I phases. Such ob- 
servations were summarized in Looney et al. (2000), who 
classified young multiple star systems into three types: (1) 
separate envelope (separation > 6000 AU), (2) common 
envelope (separation 150 — 3000 AU), and (3) common 
disk systems (separation < 100 AU). Separate envelope 
systems exhibit clearly distinct centers of gravitationally 
bound cores with separations greater than 6000 AU and 
the components are surrounded by low density material 
(e.g., NGC1333 IRS2, SVS13). Common envelope systems 
have one primary core of gravitational concentration which 
apparently breaks up into multiple fragments with sepa- 
rations of 100 - 3000 AU (e.g., IRAS 16293-2422). Com- 
mon disk systems have circumbinary disk-like envelopes 
surrounding young stellar objects with separations smaller 
than 100 AU (e.g., L1551 IRS5). 

We have also uncovered three distinct types of fragmen- 
tation of magnetic clouds in our numerical experiments. 
They correspond roughly with the above empirical classi- 



fication of young binary and multiple star systems. Based 
upon the results presented in §4, we classify theoretically 
the process of fragmentation into three types, depending 
on the initial model parameters: (1) separate core forma- 
tion, in which the cloud breaks up into two separate cores 
during the isothermal collapse phase, (2) bar fragmenta- 
tion, in which the aspect ratio of the bar becomes large 
enough to break up into multiple fragments along the ma- 
jor axis and the fragmentation takes place after the bar has 
become opaque to dust emission, and (3) disk fragmenta- 
tion, in which the central region, after becoming opaque to 
dust emission, evolves into a rapidly rotating disk, which 
fragments into multiple blobs. These cases are discussed 
in turn. 

When the initial cloud has a relatively flat central part 
and is rotating relatively rapidly (e.g., w ~ 0.2 for n ~ 4), 
it fragments into two separate cores during the isother- 
mal collapse phase in the presence of an m = 2 perturba- 
tion after the central region becomes magnetically super- 
critical (separate core formation). Morphologically, this 
case resembles the separate envelope systems discussed in 
Looney et al. (2000). The early phase of fragmentation 
also looks like the double-core in the extinction map of 
Fest 1-457 (Alves et al. 2002) (see Laundardt [2001] for 
another example of double-core, in Bok globule CB 230). 
When first produced, the dense cores have insufficient or- 
bital rotation to prevent them from spiraling towards each 
other. We expect their separation to shrink appreciably 
with time, producing perhaps a highly eccentric system. 
Whether they will eventually merge together is unclear at 
the present. We suspect that the cores will survive to be- 
come individual star/stellar system, since the density is 
much higher inside the cores than outside by the end of 
our calculations. This means that it takes the cores much 
less time to collapse onto themselves to form stars than 
to merge together. Moreover, there exists a strong mag- 
netic field outside the cores which can prevent direct core 
collision. In our separate core formation case, the den- 
sity of the central inter-core region decreases with time as 
matter falls towards the cores, whereas the field strength 
increases due to magnetic flux diffusion out of the cores 
into the inter-core region. Indeed, the central inter-core 
region may become magnetically subcritical again at the 
later stages of evolution. This region serves as a "mag- 
netic cushion" to direct core collision. We propose that 
the separate core formation is a mechanism for producing 
wide binaries with separations of ~ 10 4 AU. Each core 
will have a second chance at fragmentation at higher den- 
sities, when a first bar or rotating disk is formed. The sec- 
ondary fragmentation within each core could lead to the 
formation of a small (hierarchical) multiple star system. 
As mentioned earlier, separate core formation in the pres- 
ence of an m = 2 perturbation has much in common with 
the multiple fragmentation induced by higher order per- 
turbations that we studied previously in Li & Nakamura 
(2002). Multiple core formation is an related channel for 
wide binary and multiple star formation. 

When the initial cloud contains a mass well above the 
thermal Jeans mass and has a moderately centrally con- 
densed density distribution (e.g, r ;> 157r for n ~ 2 and 
lo ~ 0.1 — 0.2), a long bar is formed, which can break 
up into several (pairs of) fragments after the bar becomes 
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opaque to the dust emission (bar fragmentation) . The mul- 
tiple fragments grow in mass due to both merging and 
accretion. This bar fragmentation case would correspond 
to the common envelope systems of Looney et al. (2000), 
since the whole bar is embedded within a common core out 
of which the bar is formed. The bar fragmentation occurs 
at a density of ~ 10 12 cm~ 3 . According to linear per- 
turbation theory, a highly elongated bar breaks up gravi- 
tationally into fragments, with separations comparable to 
the wavelength of the most unstable linear perturbation, 
which corresponds to ~ 10 2 AU at ~ 10 12 cm~ 3 . The typ- 
ical masses of the fragments are in the substellar regime. 
Some of the fragments could be ejected out of the system 
in the course of dynamical interaction (see also Bate et al. 
2002). Those ejected may become isolated brown dwarfs 
(Reipurth & Clarke 2001) and perhaps even free-floating 
massive planets. Interestingly, the cloud parameters de- 
rived by Ciolek & Basu (2000) for the well-studied starless 
core L1544 put the cloud in the bar fragmentation cate- 
gory. Therefore, L1544 appears to be a progenitor of a 
binary or small multiple star system of the common enve- 
lope type. 

When the initial cloud mass in the central plateau region 
is comparable to the thermal Jeans mass (e.g., r <; 15n 
for n ~ 2 and u ~ 0.1), the bar mode does not grow 
much during the isothermal collapse, and it is difficult for 
the (relatively short) bar to break up into fragments (disk 
fragmentation). In this case, the central region evolves into 
a rapidly rotating disk whose size depends on the rate of 
initial rotation. The interaction of the central disk with 
the bar generates complex structures in the disk where a 
secondary bar mode is excited by rotational instability if 
(3 ;> 0.3. Typically, the size of centrifugal radius is of the 
order of 10 AU. When the disk fragments into blobs, the 
separations of the blobs should be of the same order of the 
disk radius. This case may be responsible for the forma- 
tion of binary stars with a separation smaller than 10 2 AU. 
This case would be related to the common disk system of 
Looney et al. (2000). 

As mentioned in § 2, we did not consider angular mo- 
mentum transfer in the vertical direction, which can brake 
the cloud rotation efficiently during the relatively long 
magnetically subcritical phase of evolution (Nakano & 
Tademaru 1972; Basu & Mouschovias 1994). Magnetic 
braking is thought to be responsible for the low rates of 
rotation observed in dense cores (Goodman et al. 1993) 
and adopted in this paper. Once the central region of the 
cloud becomes supercritical, the peak density increases in a 
runaway fashion, leaving little time for the magnetic brak- 
ing to operate efficiently (Basu & Mouschovias 1994). We 
therefore expect the formation of separate cores, which 
occurs in the runaway collapse phase, to be little affected 
by magnetic braking. After a bar or a disk becomes dense 
enough to be opaque to dust emission, it can be supported 



by thermal pressure for a time long compared to the local 
dynamic time, and magnetic braking can in principle be- 
come efficient again (Krasnopolsky & Konigl 2002). How- 
ever, at such high densities, magnetic fields start to de- 
couple from matter (Nishi, Nakano, & Umebayashi 1991), 
which could render the braking inefficient. We will post- 
pone a detailed treatment of the magnetic braking includ- 
ing the effects of decoupling to a future work. 

6. CONCLUSION 

In this paper, we have followed the early phases of cloud 
fragmentation in magnetic clouds, using a two-dimensional 
MHD code based on the thin-disk approximation. At the 
end of our calculations, we obtained dense, discrete frag- 
ments out of initially axisymmetric clouds perturbed by 
an m = 2 mode of modest fractional amplitude A — 0.05. 
The masses of the fragments are still small, on the order 
of 10~ 2 — 10 _1 Mq. We expect the fragments to grow to 
stellar masses, by merging and further accretion, except 
perhaps those ejected dynamically out of the system soon 
after the fragmentation. The masses of those ejected may 
remain frozen in the substellar regime. 

We found three types of fragmentation in magnetic 
clouds. They are (1) a cloud breaking up during the 
isothermal phase of cloud evolution into separate cores, 
each of which collapses more or less independently into 
a star/stellar system (panels [b] and [c] of Fig. 7), (2) 
a core collapsing into a needle- like, opaque "first bar", 
which breaks up into several fragments (Fig. 8), and 
(3) a rapidly rotating, opaque "first disk" producing self- 
gravitating blobs, which may survive to become seeds of 
relatively close binary and multiple systems (Figs. 9 and 
10). The fragments in the above cases share, respectively, 
a common cloud, common core, and common disk. They 
correspond loosely to the empirical classification of inde- 
pendent envelope, common envelope and common disk sys- 
tems for the youngest, embedded binary and multiple sys- 
tems. Much work is needed to firm up this connection. 

In the near future, we plan to extend our numerical sim- 
ulations beyond the initial fragmentation to the protostel- 
lar accretion phase, to determine the physical properties 
of the binary and multiple systems formed in magnetic 
clouds, such as their eccentricity, mass ratio, and orbital 
period. In a longer run, we hope to generalize the cal- 
culations to 3D, which will allow for a detailed study of 
magnetic braking, which may have a profound impact on 
the final orbits of the formed stellar system. 

Numerical computations in this work were carried out 
at the Yukawa Institute Computer Facilities, Kyoto Uni- 
versity. F.N. gratefully acknowledges the support of the 
JSPS Postdoctoral Fellowships for Research Abroad, and 
Z.Y.L. is supported in part by NASA. 
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Parameters of Models 
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Fig. 1 (standard model) 
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Fig. 8 (bar fragmentation) 
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Figs. 9, 10 (disk fragmentation) 
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Figs. 7b, c (separate core formation) 



Note. — All models are computed on an initial grid of 512 x 512 in the disk (x-y) plane. 
For the standard model shown in Fig. 1, the grid number is increased to 1024 x 1024 for 
each level of refinement. 
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Fig. 1. — Snapshots of the surface density distribution and velocity field for the non-rotating cloud perturbed by an m = 2 mode of fractional 
amplitude A = 0.05 at six dimensionless times: (a) t=0, (b) 18.6813, (c)25.6429, (d) 25.6849, (c) 25.6853, and (f) 25.6855. The model cloud 
has a set of "standard" parameters of (ro, n, u>, To) = (10-7T, 2, 0, 1.5). Panel (a) shows the initial state. In panels (b) through (f), only the 
central regions are shown. The contour curves in each panel are for the surface density normalized by S max , whose value is given above 
each panel. Also given are the flux-to-mass ratio (r c ) at the density peak, and length unit for each panel. The arrows are velocity vectors 
normalized by the effective isothermal sound speed c s (without magnetic contribution), whose magnitude is indicated above each panel. Note 
that in panel (a) the added perturbation lowers the flux-to-mass ratio T = 1.43 from the reference value To = 1.5 by 5%, corresponding to 
the fractional amplitude of the perturbation. For dimensional units, see § 2. 
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Fig. 2. — Evolution of (a) the aspect ratio of the highest density region with surface density above 10 1//2 E m ax, and (b) the flux-to-mass 
ratio T c at the surface density maximum. Solid curves are for the standard model with n = 2, To = 1.5, ro = 5-7T, and w = 0. 
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Fig. 3. — Distributions of surface density E, infall speed V, and the vertical magnetic field strength B z , and the two components of magnetic 
field in the disk plane (B x and B y ) along the major (solid lines) and minor (dotted lines) axes of the standard model shown in Fig. f, at 
seven different times. Note that before the breakdown of isothermal equation of state the surface density and the strength of each magnetic 
field component approach a power-law distribution of oc d^ 1 (plotted as a dashed line in panels [a], [b], and [d], with d denoting the distance 
from the center), and the infall speed is more or less constant. This self-similar behavior of isothermal collapse is explored further in §3.4. 
In panels (b) and (d), the solid and dotted lines coincide with each other at the initial state because of the initial axisymmetric magnetic 
configuration. In panel (c), the initial infall speed is equal to zero owing to the static initial condition and accordingly not depicted. 
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Fig. 4. — Snapshots of the surface density distribution and velocity field at the time when S max = 10 for the nonrotating models with 
different model parameters. Panel (a): the model with a weaker magnetic field, (ro, n, ui, To, A) = (10-7T, 2, 0, 1.25 , 0.05), showing that decreasing 
field strength promotes bar growth. Panel (b): the model with a larger ro, (ro, n,ui, To, A) = (15_7T, 2, 0, 1.5, 0.05), showing the positive effect of 
the initial cloud mass on the bar growth. Panel (c): the model with a larger n and a smaller ro, (ro, n, u, To, A) = (5tt_, 4, 0, 1.5, 0.05), showing 
the positive effect on the bar growth of a fiat density distribution. These panels should be compared with panel (d) of Fig. 1. The contours, 
arrows and notations have the same meaning as in Fig. 1. The parameters having different values from those of Fig. 1 arc underlined. 




Fig. 5. — Snapshots of the surface density distribution and velocity field of the slowly-rotating model with (ro, n, ui, To, A) = 
(107T, 2, 0.1 , 1.5. 0.05) at three different stages: (a) t = 25.025, (b) t = 25.066, and (c) t = 25.067. The contours, arrows and notations 
have the same meaning as in Fig. 1. The parameters having different values from those of Fig. 1 are underlined. Comparison with panels (c) 
through (e) of Fig. 1 points to a positive effect of the slow rotation on the bar growth. 
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Fig. 6. — Axisymmetric evolution of two clouds without perturbations. Plotted arc (a) the flux-to-mass ratio and (b) the force ratio for 
a centrally-condensed model with the standard set of cloud parameters (ro, n, ui, To) = (10-7T, 2, 0, 1.5), and (c) the flux-to-mass ratio and (d) 
the force ratio for a less centrally-condensed model with (ro, n, u>, Tq) = (5_7r, 4, 0, 1.5), at times when S c = 1.2(t = 0), 10, 10 2 , 10 3 , • • •. The 
force ratio is defined as the ratio of the effective gravity (F g e ff) to the effective pressure force (-Fp >e ff), where -F gj0 ff = -Fgrav + Ften and 
-Fp.cff = -Fgas + Fb- In these one-dimensional calculations, we do not stiffen the isothermal equation of state even when the surface density 
exceeds the critical value of E cr = 1.9 X 10 4 . In all panels, we show the distributions at the times when S c > E cr separately with dashed 
lines. 




Fig. 7. — Snapshots of the surface density distribution and velocity field of two models (a) with (ro, n, ui, To, A) = (5_7r, 4, 0.1 , 1.5, 0.05), 
and (b) and (c) with (ro, n, u, To, A) = (5_7T, 4, 0.2 , 1.5, 0.05). The contours, arrows and notations have the same meaning as in Fig. 1. The 
parameters having different values from those of Fig. 1 arc underlined. In the more rapidly rotating model, the initially-imposed m = 2 mode 
evolves into a long bar, which breaks up into two distinct cores (panel c). 
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Fig. 8. — Snapshots of the surface density distribution and velocity field of the model with (ro, n, u>, To, A) = (157T, 2, 0,2 , 1.5, 0.05) at three 
different times: (a) t = 29.5142, (b) t = 29.5144, and (c) t = 29.5145. Panels (d) though (f) show the enlargement of the central region of 
panels (a) through (c), respectively. The contours, arrows and notations have the same meaning as in Fig. 1. The parameters having different 
values from those of Fig. 1 are underlined. In this model, the aspect ratio of the first bar becomes very large and the fragmentation starts 
near the center and at the ends (panel a). Subsequent fragmentation occurs successively in the high density central region, and seemingly 
propagates towards the ends (panel b). At the same time, the central fragments merge together to form more massive fragments (panel b). 
This process is repeated at the time shown in panel (c). 
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Fig. 9. — Snapshots of the surface density distribution and velocity field of the model with (ro, n, u>, To, A) = (5n, 2. 0.1, 1.5, 0.05) at six 
different times. The contours, arrows and notations have the same meaning as in Fig. 1. The parameters having different values from those 
of Fig. 1 are underlined. As the infall along the major axis becomes significant, the truly first core is formed at the center, which evolves into 
a rapidly rotating disk (panel a). In the rotating disk, an m = 2 mode develops, which interacts with the inflows along the first bar (panels 
b and c), generating shock waves in several places (panels d and e). The disk fragments into several blobs due to complex shock interaction 
(panels c and f). 
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Fig. 10. — Surface density distributions of the central rotating disk for the model with (ro, n, u>, Tq, A) = (57T, 2, 0.1. 1.5. 0.05) at three 
different times. Panels (a), (b) and (c) are at the same times as panels (b), (d), and (f) of Fig. 9. The surface density contours are linear in 
scale, plotted at levels of 0.1, 0.2, . . ., 0.9 of the maximum, Lmax' At the time shown in panel (c), the bar breaks up into two blobs, which 
appear to be self-gravitating. 



